Excitonic effects in solids described by time-dependent density functional theory 
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Starting from the many-body Bethe-Salpeter equation we derive an exchange-correlation kernel 
fxc that reproduces excitonic effects in bulk materials within time-dependent density functional 
theory. The resulting fxc accounts for both self-energy corrections and the electron-hole interaction. 
It is static, non-local and has a long-range Coulomb tail. Taking the example of bulk silicon, we show 
that the —a/q^ divergency is crucial and can, in the case of continuum excitons, even be sufficient 
for reproducing the excitonic effects and yielding excellent agreement between the calculated and 
the experimental absorption spectrum. 
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The calculation of electronic excitations has remained 
a major challenge in the field of solid state theory. In fact, 
whereas ground state properties can be computed to- 
day with good precision within Density Functional The- 
ory (DFT) [0, electronic excitations are accessible only 
through additional corrections. Within Many-Body Per- 
turbation Theory, Hedin's GW corrections |^ are used 
to get electron addition and removal energies, and the 
Bethe-Salpeter Equation (BSE) for neutral excitations 
as those measured for example in absorption or elec- 
tron energy loss spectroscopy (EELS). In fact, since sev- 
eral years ab initio BSE calculations yield generally good 
agreement between the calculated and the experimental 
absorption spectra for both finite and infinite systems 
1^-^ . However these calculations are necessarily cumber- 
some, which has up to now prevented their application 
to very complex systems. 

In the last years, an alternative approach has been 
developed which allows in principle to correctly de- 
scribe exchange-correlation effects in the neutral excited 
state, namely Time-Dependent Density Functional The- 
ory (TDDFT) [01]. As in the case of static DFT the 
main obstacle resides in finding a good approximation 
to the unknown exchange-correlation (xc) contribution. 
For the ground state properties of the majority of fi- 
nite and infinite systems, the local density approximation 
(LDA) has turned out to yield surprisingly good results. 
In the case of the absorption and EELS spectra of fi- 
nite systems, and of EELS spectra of solids, TDDFT has 
yielded good results using the adiabatic LDA approxima- 
tion (TDLDA) for the xc-kernel fxc However, this is 
not true for the absorption spectra of solids In fact, 
it has turned out that the results for the latter obtained 
within TDLDA are extremely close to those obtained in 
a simple Random Phase Approximation (RPA) calcula- 
tion, where the xc-kernel is completely neglected, and 
only local field effects (in other words, the contribution 
coming from the variation of the Hartree potential) are 
taken into account. It should be pointed out that this 
similarity between RPA and TDLDA spectra to some ex- 
tent also holds in the case of clusters and EELS spectra 
of solids, but the results are nevertheless satisfactory be- 
cause in those cases already RPA spectra, including local 
field effects, are in good agreement with experiment. 

On the other hand, it would be extremely desirable to 
obtain good absorption spectra of solids within TDDFT, 
since the equations to be solved are two-point ones, in 
contrast to the four-point BSE. The obstacle to be re- 
moved to this aim is hence the fact that a good approxi- 
mation to the xc-kernel must be found. In this context, 
extensive discussions can be found in literature about the 
need to include long-range nonlocal terms and dynamical 
(memory) effects in the kernel PJTo|. 

In this work, we show how a TDDFT equation for the 
macroscopic dielectric function can be derived from the 
Bethe-Salpeter equation. The derivation is exact in the 
sense that the resulting equation does yield the same 
spectra as the standard Bethe-Salpeter equation (with 



its own various approximations [^|-^). We demonstrate 
that the resulting xc-kernel has a contribution, the 
strength of which is inversely proportional to the screen- 
ing in the system. Finally, we calculate the absorption 
and the refraction index spectra of bulk silicon within 
TDDFT, using only this static long range tail as xc- 
kernel, and obtain excellent agreement with experiment. 

To start with it is useful to put the TDDFT equa- 
tion and the Bethe-Salpeter equation for e"^ on the same 
footing. In fact, both equations can be written schemat- 
ically in the same Dyson-like form. 



S = ^(0) + S^°^KS. 



(1) 



Here, S can be either the two-point polarisability Xj from 
which we can obtain the inverse dielectric function e^^ = 
1 -|- vx, or S = L, the two-particle correlation function 
which yields x (and then e~^) by contracting two of its 
four indices 



x{xi,X2) = L{xi,xf;x2,x^). 



(2) 



Here x stands for space, spin and time coordinates. Of 
course this holds only for what concerns the form, but not 
the specific details. First, quantities in the TDDFT equa- 
tion are two-point ones whereas in the BSE they are four- 
point ones. Second, in TDDFT 5'^°^ is the independent- 
particle response function x'"' constructed with the 
Kohn Sham (KS) orbitals and eigenvalues, whereas in the 
BSE formalism S'^"^ stands for the independent quasipar- 
ticle response , i.e., it is constructed using quasiparti- 
cle eigenvalues and eigenf unctions (e.g., obtained from a 
GW calculation). Finally, in the TDDFT case, the ker- 
nel K is defined as K = v + ^tddft^ where v is the 
Coulomb potential and ^^tddft s|;a,nds for the fxc ker- 
nel. In the case of the BSE, K = v + F^^^, where v 
and i^BSE g^j^g ^Yie four-point functions v{xi,X2,X3,X4) = 
5{xi,X2)5{x3,X4)v{xi,X3) and F^^^{xi, X2, X3, x^) = 
-~S{x2j X4)S{xi,X3)W{xi, X2), where W is the screened 
interaction. In the case of TDDFT, one can also un- 
derstand Eq. (|l|) as a four-point equation, but with 

^TDDFT^^^^2;2,X3,X4) = 6{xi, X2)S{x3, X4) fxc{xi, X3) 

which implies that one can immediately contract the in- 
dices by pairs and reduce the equation to a two-point 
one. 

If one performs a basis transformation of the form 
xi — * ipnixi), where the ipnixi) are the one-particle or- 
bitals which diagonalise the four-point 5''-*'^ , Eq. (|l|) can 
be transformed to an effective two-particle Hamiltonian 
equation, 



H 



{nin2){n3ni) — (^"2 ) <^niri3 ^Ji2n4 



~^ (/"l /n2 ) -^(nin2)(n3"4) ■ 



(3) 



Here the indices Ui refer to the fact that matrix elements 
involve four eigenfunctions of the starting effective one- 
particle Hamiltonian with eigenvalue e„. and occupation 
fm- Defining the identity operator / = Smim3Sm2m4,, S 
is directly obtained from H^^ as 
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'(nin2)(n3"4) 



)(n3n4 



ifn.-fns)- (4) 



can be diagonalized, and from its eigenvalues E\ and 
eigenstates A^^^^ the spectral representation of S can be 
constructed. 

Since the ipn{xi) have to diagonalise S'^'^\ they must 
be the KS orbitals for the TDDFT equation and the 
QP eigenfunctions for the BSE. If these functions are 
equal, and if the and Ex are equal, the BSE and 
TDDFT spectra would be the same. For a static f^c ker- 
nel, this implies that the matrix elements of the Hamilto- 
nians -ffxDDFT ^'^^ ^bse equal. In this scenario, we 
can directly compare the BSE and TDDFT approaches. 
First, in the BSE the eigenvalues e„ are quasiparticle en- 
ergies (as obtained for example from a GW calculation), 
whereas in TDDFT they are the eigenvalues obtained 
from the KS equation with the (in principle exact) xc- 
potential. Second, for the BSE we have 



F, 



BSE 



j c?r(ir'$(ni,n3;r)H/(r',r)cf>*(n2,n4;r') 



with the matrices $ defined according to 

$(ni,n2;r) := V„i(r)V':,(r), (5) 

whereas the exchange-correlation contribution to the 
TDDFT kernel reads 



F, 



TDDFT 

ni"2)(n3n4) 



= j c^rdr'$(nl,n2;r)^c(r,r')$*(^^3,?^4;r') 



(The Hartree contribution 'f(„j„2)(„3„^) is of course equal 
in both cases) . A comparison of the two cases does hence 
immediately tell us that the BSE and TDDFT equations 
would yield the same spectrum if the static fxc satisfies 

(/«i-/«2) j drdr'$(ni,n2;r)/:rc(r,r')$*(n3,n4;r') = 

■F(nin2)(n3n4,)j (6) 

where 

T — /'cQP _ eQP _ e^FT , DFT\ r r 

(nin2)(n3ri4) — \'^n2 ^ni ^n2 ' ^rii ) ''nin3^n2ni 

+ (/ni ~ /n2)-P'(nin2)(n3Ti4)- C^) 

It is clear that, if the transformation xi — > ipn{xi) was 
complete in all four indices, Eq. (^) could never be sat- 
isfied. The reason is that the two operators, J^tddft 
and cannot he equal because of the way the 5 func- 
tions are put in real space. On the other hand, if the 
two operators cannot be made equal, then the spectrum 
can only be equal if at least one of the two operators (in 
that case, fxc) is energy dependent [ |12| . This in principle 
correct, general statement, can be made less restrictive 
by realizing that in practice only a finite number of tran- 
sitions contributes to the optical spectrum. This means 
that we can use an incomplete basis in transition space. 
In this reduced Hilbert space we can still find a static 



operator that satisfies the required equality in transition 
space in a particular energy range, even though the real- 
space operators are not equal In order to discuss 
this possibility, we rewrite Eq. ^ as 



(/«! -/«2) E 1'("i,"2;G)/,c(q,G,G')<i>*(n3,n4;G') 



G,G' 



•^(nin2)(ri3Ti4) J 



(8) 



where q = k2 — ki = k4 — k3. Since for the particle- hole 
and hole-particle contributions of a non-metal the factor 
(/«i ~ fn2) can never be zero, we could in principle obtain 
the matrix f^c from Eq. (p|) by inverting the matrices 



/,,(q,G,G')= 



1 



(/ni fn2) 



$ ^(ni,n2;G)jr(„^„^)(„3„^)($*) 1(^3, 714; G'). (9) 

As pointed out above, the sum over the indices is 
necessarily limited to some subspace of important transi- 
tions, because otherwise the inverse of $ does not exist. 
In certain cases, a subspace which at the same time allows 
the matrix to be invertible and to reasonably reproduce 
the exciton spectrum cannot be found. One example is 
a two-band model with g^FT _ ^qp consisting of plane 
wave states ipbki^) = e^^'''^ e^^'^ , with Gf, = G^, or Gc- In 
this case the matrices $ are k-independent, which means 
that all their rows are equal. In other words, the ma- 
trix i^TDDFT k-independent, whereas the matrix i^^SE 
goes as l/(fc — fc')^, which is a clear contradiction. On 
the other hand, if a subspace where $ is invertible can 
be found, it is clear that the resulting kernel is frequency 
independent, as none of the quantities implied, J'^S^ or 
$, are frequency dependent (if the kernel of the BSE is 
chosen to be frequency independent fl^). Moreover f^c 
is necessarely not local, as the expression in Fourier space 
depends separately on G and on G'. In fact, its nonlo- 
cality is used as the degree of freedom which is necessary 
to fulfill equality (^. 

On the other hand, one can try to look directly at S 
instead of focusing on f^c- To do that we go back to the 
initial Eq.(^ which we can write in a symmetric way as 



5^5(o)(^(o)_^(o)^^(o))-i^(o)^ 



(10) 



Since the independent particle polarisability is 



gDFT _ gDFT _ ^ 



the equation for the two-point S — x becomes 



X(r,r',w) 



dr"dr'"x^°\v. 



X^°\ri,r2, uj) 



rfi"3dr4 x^°' (ri , r3, U})v{r3,r4)x^°^ (r4, r2, w) 
^ $*(ni,n2;ri) 



,DFT 



,DFT 
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^ ^(n3,n4;r2)(/n4 - fnj ] 

(nin2)(n.3"4) ,DFT „ ,DFT _ , , 

'^n4 "=713 ^ Jr",r"' 

■X^'\r"'y,u;), (11) 

where we have used ^{ni, n2, r) — $*(n2, rti, r), and sub- 
stituted the term ^fxc^* using Eq. (^. In other words, 
for those cases where the equality ^ can be fulfilled, i.e. 
in particular when a static kernel can be found, we have 
succeeded in writing the two-point TDDFT equation in 
a way that exactly yields the BSE spectrum Q. 

The explicit knowledge of the kernel is actually not 
needed. There is still a four-point quantity appearing, 
namely however only in a matrix product instead of 
inversions and diagonalisations, which allows to change 
space in a convenient way as it is also done in recur- 
sive inversions of BSE 1). Despite this fact, in view of 
the ongoing discussions about the xc-kernel it is inter- 
esting to examine some of its features, and in particular 
its long-range behaviour. In fact, for valence (v,k) and 
conduction (c,k-|-q) states, $(?;, k, c, k + q; G = 0) goes 
to zero as q for small q. Since J'{v.c).{v,c) in this limit 
behaves as a constant, an fxc{l, G = G' = 0) obtained 
from Eq. must behave as l/q^. There is in fact a 
positive long-range contribution stemming from the QP 
shift of eigenvalues (as also predicted in ref. [|l5|), and a 
negative one resulting from the electron-hole interaction, 
which is the main point of interest here. 

Using the above results, one can also understand why 
the TDLDA approximation yields much worse results for 
the absorption than for the loss spectra of solids. In fact 
energy-loss spectra are directly related to the inverse of 
the dielectric matrix = 1-1- For the calculation 
of X; the kernel fxc is then added to u, which already 
contains a long-range contribution v{G = 0). So the 
presence or absence of the long range term in fxc does not 
necessarily show up. However, in the case of absorption 
spectra, one can show |l^ that the macroscopic dielectric 
function is given by 

eM(t^) = 1 - lim[i;(g)xG=G'=o(q,t^)], (12) 

q-»0 

where x has been calculated in a Dyson-like equation 
such as using the same fxc, but added to a coulombian 
V which does not contain the long range term, i.e. v{G = 
0) is set to zero. Obviously in that case, a neglect of the 
divergence in fxc makes an essential difference. 

We can carry this discussion about the long-range 
electron-hole interaction term further by a) assuming 
that we absorb the first, positive contribution in the 
energy shift of our starting (since anyway we do 
not know the eigenvalues of the exact xc-potential which 
would go along with the exact kernel) and b) sup- 
posing that we have a system where the long-range 
term is completely dominating the rest of the xc con- 
tribution, namely, where we can approximately write 
fxc{<l,G,G') — — (5G,G'a/|q + Gp. In other words, the 
long-range electron-hole attraction part of fxc reduces 



the Hartree part, which is reasonable since it stems from 
the exchange potential. This approximation works best 
for systems with weakly bound excitons. For a demon- 
stration we have therefore performed a TDDFT calcu- 
lation for bulk silicon, in the following way: first, we 
have determined the DFT-LDA electronic structure. Sec- 
ond, we have constructed x^^\ but with the eigenvalues 
shifted to the GW ones, in order to simulate the first 
part of the kernel as explained above. Third, we have 
used /a;c(r, r') = — a/47i"|r — r'|, with the empirical value 
a = 0.2. The result of the TDDFT calculation for euito) 
is shown in Fig. 1. The dots are the experimental results 
for the absorption spectrum (Im(£)) measured by Laut- 
enschlager et al. [|l^ and the refraction index (Re(e)) 
measured by Aspnes and Studna jlj]. The dot-dashed 
curve is the result of a standard TDLDA calculation (i.e. 
using DFT-LDA eigenvalues and the static short-range 
LDA xc kernel). We find the well-known discrepancies 
with experiment. The dashed curve is the BSE result. 
Finally, the continuous curve is the result of our approx- 
imate TDDFT calculation: it fits almost perfectly all 
experimental features in both real and imaginary parts 
of s. It turns hence out that this static long-range contri- 
bution to the kernel is sufficient to reproduce the strong 
excitonic effect in a material with weakly bound excitons 
like silicon. 

In conclusion, we have derived a TDDFT equation 
from the Bethe-Salpeter equation which should be partic- 
ularly suitable for practical applications to the absorption 
spectra of solids. We have demonstrated that the static 
exchange-correlation kernel has a long-range contribution 
stemming from the electron-hole interaction. We have ex- 
plained why this long-range contribution is particularly 
important for the absorption spectra of solids. At the ex- 
ample of bulk silicon, we have shown how a very simple 
approximation for the kernel can yield excellent agree- 
ment between the calculated TDDFT absorption spec- 
trum and experiment. 
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FIG. 1. Silicon, optical absorption (bottom) and refraction index (top panel) spectra. Dots: experiment. Dot-dashed curve: 
TDLDA result. Dashed curve: result obtained through the Bethe-Salpeter method. Continuous curve: TDDFT result using 
the long-range kernel derived in this work. 
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